TIRE-seq Dendritic cells Flt3 stimulation

Recap

Prior 96w evaluation of TurboCapture-Seq v2 showed low UMIs recovered and high seq saturation. I did low throughput troubleshooting and didn’t see any issues. Process this experiment myself taking care to remove residual liquids from wash steps.

Process Hui Shi of Naik lab Bcor + Flt3 timecourse. Includes a few of my samples.

Samples

Sorted dendritic cells
HEK293T cell lysates in 1x TCL @ cells/uL
PBMC cell lysates in 1x TCL @ 500 cells/uL
No template control 1x TCL

Notebook recap

SCE object in generate in 1A_generateSCE_reads notebook. The samples were then clustered in the 2B Clustering Wt notebook

Notebook aim

Compare the mixed cell culture stimulated with Flt3 for one day. This is a basic contrast that should have a lot of differentially expressed genes.

Workflow is from:
https://bioconductor.org/packages/devel/workflows/vignettes/RNAseq123/inst/doc/limmaWorkflow.html#differential-expression-analysis

Read SCE and preprocessing

This was generated in notebook 2B.

sce <- readRDS(here::here(
   "data/TIRE_dendritic_mouse/SCEs/DCs_cluster.sce.rds"))

sce <- sce[,sce$Cell_type == "Total_DC_culture"]

dge <- scran::convertTo(sce, type="edgeR")

Have a look at the important metadata.
There are not enough replicates for the preDC contrast to be meaningful.

tb <- dge$samples[,c("Ligand", "Timepoint_Day")]

tb %>% 
  count(Ligand, Timepoint_Day)

Recap the PCA

Samples are distinct but group along PC2

plt2 <- scater::plotPCA(sce,colour_by="Ligand") + 
  theme_Publication()

plt2

The samples are treated with FLt3 ligand and harvest after 24 hours.

table(dge$samples[,c("Timepoint_Day", "Ligand")])

##              Ligand
## Timepoint_Day Flt3L None
##             0     0    6
##             1     3    0
##             5     0    0
##             7     0    0

Filter low expressed genes

I added this last to filter genes. About half the genes are filtered out here.
Doing this reduces the multiple testing burden and fits variation better.

dim(dge)

## [1] 21484     9

keep.exprs <- filterByExpr(dge, group=dge$samples$Timepoint_Day, min.count=1)
dge <- dge[keep.exprs,, keep.lib.sizes=FALSE]
dim(dge)

## [1] 9550    9

Design matrix

Build the model matrix.
Regress out the effect of the timepoint as this is minor compared to the cell type subsets.

This means I can explicitly define a contrast matrix to make the comparisons of interest.

sm <- model.matrix(~0+Ligand, data=dge$samples)
head(sm)

##            LigandFlt3L LigandNone
## ACAGTAGCTC           0          1
## AGACGCATCG           1          0
## CCGTATGTAG           0          1
## CGTCAATAGT           0          1
## CTTGATCGCG           0          1
## GGAAAGATAC           0          1

# hypens not allowed
colnames(sm) <- make.names(colnames(sm), unique = FALSE, allow_ = TRUE)

Decide on the contrasts. In this case its simply plus or minus Flt3 ligand.

contr.matrix <- makeContrasts(
   Ligand = LigandFlt3L - LigandNone, 
   levels = colnames(sm))

contr.matrix %>% 
  kable()

	Ligand
LigandFlt3L	1
LigandNone	-1

In each contrast, the format is A - B where:

A represents the condition considered as the “treatment” or point of interest
B represents the condition considered as the “control” or baseline

Limma voom

Fit this model using limma voom

par(mfrow=c(1,2))
v <- voom(dge, sm, plot=TRUE)

vfit <- lmFit(v, sm)
vfit <- contrasts.fit(vfit, contrasts=contr.matrix)
efit <- eBayes(vfit)
plotSA(efit, main="Final model: Mean-variance trend")

Differential expression analysis

Check how many genes are differentially expressed

summary(decideTests(efit))

##        Ligand
## Down      313
## NotSig   8717
## Up        520

Check what the genes are

ligand <- topTable(efit, coef=1, n=length(efit$genes$ID), sort.by = "logFC")
head(ligand)

ChatGPT explanations:

JunB in particular is a key transcriptional regulator of dendritic cell maturation and activation.
JunB was found to bind extensively to the chromatin of LPS-treated dendritic cells and regulate the expression of genes involved in the inflammatory response.
Chad encodes a cartilage matrix protein called chondroadherin, which is thought to mediate adhesion of isolated chondrocytes (cartilage cells).
- It promotes attachment of chondrocytes, fibroblasts, and osteoblasts (bone cells)
- SO CHad doesn’t make aa alot of sense but this is a mixed cell culture
FosB and its truncated isoform delta-FosB have also been implicated in the maturation of human monocyte-derived dendritic cells
Zfp36 (Zinc finger protein 36):
- An RNA-binding protein that plays a crucial role in regulating inflammation and immune responses. It binds to AU-rich elements in the 3’ untranslated regions (UTRs) of target mRNAs, leading to their degradation.
- This post-transcriptional regulation helps control the production of pro-inflammatory cytokines such as TNF-alpha.
- Zfp36 helps modulate the inflammatory response and maintain immune homeostasis by controlling the levels of cytokines and other inflammatory mediators.
Myl10 (Myosin Light Chain 10):
- Regulatory light chain of myosin, which is involved in muscle contraction and cell motility.
- Important for cytoskeletal dynamics, which are crucial for cell migration, antigen uptake, and the formation of immune synapses. Proper regulation of the cytoskeleton enables dendritic cells to efficiently travel to lymph nodes and present antigens to T cells.
Flt1 (Fms-Related Receptor Tyrosine Kinase 1):
- Also known as VEGFR-1 (Vascular Endothelial Growth Factor Receptor 1), is a receptor tyrosine kinase that binds to VEGF (Vascular Endothelial Growth Factor).
- Flt1 plays a role in the migration and function of dendritic cells by influencing the vascular environment and facilitating their movement through tissues.

Visualise differential expression testing

results <- as_tibble(ligand)
results$ID <- rownames(ligand)

# add a column of NAs
results$DElabel <- "n/s"
# if log2Foldchange > 1 or < -1 and pvalue < 0.05, set as "UP"
results$DElabel[results$logFC > 1 & results$adj.P.Val < 0.1] <- "Flt3"
results$DElabel[results$logFC < -1 & results$adj.P.Val < 0.1] <- "Ctrl"

Add gene labels

results$genelabels <- ""
results$genelabels[results$ID == "Fosb"] <- "Fosb"
results$genelabels[results$ID == "Jun"] <- "Jun"
results$genelabels[results$ID == "Chad"] <- "Chad"
results$genelabels[results$ID == "Zfp36"] <- "Zfp36"
results$genelabels[results$ID == "Egr1"] <- "Egr1"
results$genelabels[results$ID == "Ier2"] <- "Ier2"

results$genelabels[results$ID == "Fos"] <- "Fos"
results$genelabels[results$ID == "Nfkbiz"] <- "Nfkbiz"
results$genelabels[results$ID == "Ftl1"] <- "Ftl1"
results$genelabels[results$ID == "Lpl"] <- "Lpl"
results$genelabels[results$ID == "Siglech"] <- "Siglech"
results$genelabels[results$ID == "Myl10"] <- "Myl10"

Volcano

library(ggrepel)

plt3 <- ggplot(data=results, aes(x=logFC, y=-log10(adj.P.Val), colour=DElabel, label=genelabels)) + 
  geom_point(alpha=0.33, size=1.5) +
  geom_text(size=3.5, alpha=1, colour="black",nudge_y=0.5) +
  guides(colour = guide_legend(override.aes = list(size=3, alpha=1))) +
  scale_color_manual(values = c("darkblue", "grey", "darkorange"), name = "Upregulated") +
  geom_vline(xintercept = 1, linetype="dotted") + geom_vline(xintercept = -1, linetype="dotted") +
  theme_Publication()

plt3

MAplot

plt2 <- ggplot(data=results, aes(x=AveExpr, y=logFC, colour=DElabel)) + 
  geom_point(alpha=0.75, size=1.5) +
  guides(colour = guide_legend(override.aes = list(size=3, alpha=1))) +
  ylab("Log fold chnage") + xlab("Log counts per million") +
  scale_color_manual(values = c("darkblue", "grey", "darkorange"), name = "Upregulated") +
  theme_Publication()

plt2

Gene set testing with camera

Need to convert geneIDs from ensembl to enterez

geneids <- as.data.frame(v$genes$ID)
colnames(geneids) <- "ENSEMBL"

geneids$entrez <- mapIds(org.Mm.eg.db, keys = geneids$ENSEMBL, keytype = "ENSEMBL", column = "ENTREZID")

Look at C2 curated genesets

load("data/MSigDB/mouse_c2_v5p2.rdata")
idx <- ids2indices(Mm.c2,identifiers = geneids$entrez)
cam.Ligand <- camera(v,idx,sm,contrast=contr.matrix[,1])
head(cam.Ligand,10)

par(mfrow=c(1,1))

barcodeplot(efit$t[,1], index=idx$REACTOME_CHOLESTEROL_BIOSYNTHESIS, 
            index2=idx$NAGASHIMA_EGF_SIGNALING_UP, main="Ligand")

(top) NAGASHIMA_EGF_SIGNALING_UP (bottom) REACTOME_CHOLESTEROL_BIOSYNTHESIS

Generate gene set barchart

geom_GeneSet_Barchart(cam.Ligand)

Try H hallmark gene sets

load("data/MSigDB/mouse_H_v5p2.rdata")
idx <- ids2indices(Mm.H,identifiers = geneids$entrez)
cam.Ligand <- camera(v,idx,sm,contrast=contr.matrix[,1])
head(cam.Ligand,10)

The Hallmark genesets are more informative to me.

par(mfrow=c(1,1))

barcodeplot(efit$t[,1], index=idx$HALLMARK_CHOLESTEROL_HOMEOSTASIS, 
            index2=idx$HALLMARK_TNFA_SIGNALING_VIA_NFKB, main="Ligand")

(top)HALLMARK_TNFA_SIGNALING_VIA_NFKB (bottom)HALLMARK_MTORC1_SIGNALING